#include "ProgressiveQD.h"
#include "Utils.h"

ProgressiveQD::ProgressiveQD( std::shared_ptr<TaylorSeries> taylorSeries, double c0 ) : _taylorSeries(taylorSeries), _c0(c0), _n(1)
{
}

ProgressiveQD::~ProgressiveQD(void)
{
}

double ProgressiveQD::d( unsigned int n )
{
	if(n < _vd.size())
		return _vd[n];

	if(n == 0){
		_vd.push_back(1/c(0));
		return _vd[0];
	}

	double sum = 0.0;
	for(unsigned int k = 1; k <= n; k++)
	{
		sum += c(k)*d(n-k);
	}

	_vd.resize(n+1, 0.0);
	_vd[n] = (-1/c(0))*sum;

	return _vd[n];
}

double ProgressiveQD::c( unsigned int n )
{
	double coeff = _taylorSeries->coefficient(n); 
	
	// the algorithm requires all coefficients are non zero
	if( isEqual(coeff,0) ){ 
		if(n==0)
			return _c0; // special handling for nonunit series
		else
			throw std::logic_error("All coefficients must be non-zero in order to apply the progressive Q.D. algorithm");
	}

	return coeff;
}

double ProgressiveQD::e( int k, int l )
{
	auto it = _me.find(k);

	if(it == _me.end())
		it = _me.insert(pair<size_t, vector<double> >(k, vector<double>())).first;

	auto& vek = it->second;
	if(l < vek.size())
		return vek[l];

	double val = (q(k-1,l+1)*e(k-1, l))/q(k,l);

	vek.push_back(val);

	return val;
}

double ProgressiveQD::q( int k, int l )
{
	auto it = _mq.find(k);

	if(it == _mq.end())
		it = _mq.insert(pair<size_t, vector<double> >(k, vector<double>())).first;

	auto& vqk = it->second;
	if(l < vqk.size())
		return vqk[l];

	double val = (e(k-1,l) - e(k, l-1)) + q(k-1,l);

	vqk.push_back(val);

	return val;
}

shared_ptr<vector<double> > ProgressiveQD::partialNumerators( int n )
{
	auto numerators = make_shared<vector<double> >();
	if(numerators->size() == 0)
		numerators->push_back(c(1));

	initialize(n);

	// iterate through the diagonals
	for(int i = _n-2; i <= n-2; i++)
	{
		double parNum = 0.0;
		// iterate through the entries in the diagonal
		for(int j = 0; j <= i; j++)
		{
			int k = -(n-2)+i;
			int l = (n-1) - i + floor(double(j)/2);

			if(j % 2 == 0)
				parNum = q(k+1, l);	
			else
				parNum = e(k+1,l);

			if(k == 0)
				numerators->push_back(-parNum);
		}
	}

	return numerators;
}

void ProgressiveQD::initialize(int n)
{
	// note: the inserts will fail if the entries have already been added
	double q0[2] = {0.0, -(d(1)/d(0))};
	_mq.insert(make_pair(0, vector<double>(q0, q0 + 2)));
	_me.insert(make_pair(1, vector<double>(1, 0.0) ));
	double e0[2] = {0.0, d(2)/d(1)};
	_me.insert(make_pair(0, vector<double>(e0, e0 + 2)));
	
	for(int k = _n; k <= n; k++)
	{
		vector<double> vq(k+2,0.0);
		_mq.insert(pair<size_t, vector<double> >(-k, vq));

		vector<double> ve(k+2,0.0);
		ve[k+1] = d(k+2)/d(k+1);
		_me.insert(pair<size_t, vector<double> >(-k, ve));
	}
	_n = n;
}
